/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkGradientFilter.cxx

  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notice for more information.

=========================================================================*/

/*----------------------------------------------------------------------------
 Copyright (c) Sandia Corporation
 See Copyright.txt or http://www.paraview.org/HTML/Copyright.html for details.
----------------------------------------------------------------------------*/

#include "vtkGradientFilter.h"

#include "vtkCell.h"
#include "vtkCellData.h"
#include "vtkCellDataToPointData.h"
#include "vtkDataArray.h"
#include "vtkDataSet.h"
#include "vtkIdList.h"
#include "vtkImageData.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkMath.h"
#include "vtkNew.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkPolyData.h"
#include "vtkRectilinearGrid.h"
#include "vtkSmartPointer.h"
#include "vtkStreamingDemandDrivenPipeline.h"
#include "vtkStructuredGrid.h"
#include "vtkUnstructuredGrid.h"

#include <limits>
#include <vector>

//------------------------------------------------------------------------------

vtkStandardNewMacro(vtkGradientFilter);

namespace
{
// special template macro for only float and double types since we will never
// have other data types for output arrays and this helps with reducing template expansion
#define vtkFloatingPointTemplateMacro(call)                                                        \
  vtkTemplateMacroCase(VTK_DOUBLE, double, call);                                                  \
  vtkTemplateMacroCase(VTK_FLOAT, float, call)

// helper function to replace the gradient of a vector
// with the vorticity/curl of that vector
//------------------------------------------------------------------------------
template <class data_type>
void ComputeVorticityFromGradient(data_type* gradients, data_type* vorticity)
{
  vorticity[0] = gradients[7] - gradients[5];
  vorticity[1] = gradients[2] - gradients[6];
  vorticity[2] = gradients[3] - gradients[1];
}

template <class data_type>
void ComputeDivergenceFromGradient(data_type* gradients, data_type* divergence)
{
  divergence[0] = gradients[0] + gradients[4] + gradients[8];
}

template <class data_type>
void ComputeQCriterionFromGradient(data_type* gradients, data_type* qCriterion)
{
  // see http://public.kitware.com/pipermail/paraview/2015-May/034233.html for
  // paper citation and formula on Q-criterion.
  qCriterion[0] =
    -(gradients[0] * gradients[0] + gradients[4] * gradients[4] + gradients[8] * gradients[8]) /
      2. -
    (gradients[1] * gradients[3] + gradients[2] * gradients[6] + gradients[5] * gradients[7]);
}

// Functions for unstructured grids and polydatas
template <class data_type>
void ComputePointGradientsUG(vtkDataSet* structure, vtkDataArray* array, data_type* gradients,
  int numberOfInputComponents, data_type* vorticity, data_type* qCriterion, data_type* divergence,
  int highestCellDimension, int contributingCellOption);

int GetCellParametricData(
  vtkIdType pointId, double pointCoord[3], vtkCell* cell, int& subId, double parametricCoord[3]);

template <class data_type>
void ComputeCellGradientsUG(vtkDataSet* structure, vtkDataArray* array, data_type* gradients,
  int numberOfInputComponents, data_type* vorticity, data_type* qCriterion, data_type* divergence);

// Functions for image data and structured grids
template <class Grid, class data_type>
void ComputeGradientsSG(Grid output, vtkDataArray* array, data_type* gradients,
  int numberOfInputComponents, int fieldAssociation, data_type* vorticity, data_type* qCriterion,
  data_type* divergence);

bool vtkGradientFilterHasArray(vtkFieldData* fieldData, vtkDataArray* array)
{
  int numarrays = fieldData->GetNumberOfArrays();
  for (int i = 0; i < numarrays; i++)
  {
    if (array == fieldData->GetArray(i))
    {
      return true;
    }
  }
  return false;
}

// generic way to get the coordinate for either a cell (using
// the parametric center) or a point
void GetGridEntityCoordinate(
  vtkDataSet* grid, int fieldAssociation, vtkIdType index, double coords[3])
{
  if (fieldAssociation == vtkDataObject::FIELD_ASSOCIATION_POINTS)
  {
    grid->GetPoint(index, coords);
  }
  else
  {
    vtkCell* cell = grid->GetCell(index);
    double pcoords[3];
    int subId = cell->GetParametricCenter(pcoords);
    std::vector<double> weights(cell->GetNumberOfPoints() + 1);
    cell->EvaluateLocation(subId, pcoords, coords, &weights[0]);
  }
}

template <class data_type>
void Fill(vtkDataArray* array, data_type vtkNotUsed(data), int replacementValueOption)
{
  switch (replacementValueOption)
  {
    case vtkGradientFilter::Zero:
      array->Fill(0);
      return;
    case vtkGradientFilter::NaN:
      array->Fill(vtkMath::Nan());
      return;
    case vtkGradientFilter::DataTypeMin:
      array->Fill(std::numeric_limits<data_type>::min());
      return;
    case vtkGradientFilter::DataTypeMax:
      array->Fill(std::numeric_limits<data_type>::max());
      return;
  }
}
template <class data_type>
int GetOutputDataType(data_type vtkNotUsed(data))
{
  if (sizeof(data_type) >= 8)
  {
    return VTK_DOUBLE;
  }
  return VTK_FLOAT;
}
} // end anonymous namespace

//------------------------------------------------------------------------------
vtkGradientFilter::vtkGradientFilter()
{
  this->ResultArrayName = nullptr;
  this->DivergenceArrayName = nullptr;
  this->VorticityArrayName = nullptr;
  this->QCriterionArrayName = nullptr;
  this->FasterApproximation = 0;
  this->ComputeGradient = 1;
  this->ComputeDivergence = 0;
  this->ComputeVorticity = 0;
  this->ComputeQCriterion = 0;
  this->ContributingCellOption = vtkGradientFilter::All;
  this->ReplacementValueOption = vtkGradientFilter::Zero;
  this->SetInputScalars(
    vtkDataObject::FIELD_ASSOCIATION_POINTS_THEN_CELLS, vtkDataSetAttributes::SCALARS);
}

//------------------------------------------------------------------------------
vtkGradientFilter::~vtkGradientFilter()
{
  this->SetResultArrayName(nullptr);
  this->SetDivergenceArrayName(nullptr);
  this->SetVorticityArrayName(nullptr);
  this->SetQCriterionArrayName(nullptr);
}

//------------------------------------------------------------------------------
void vtkGradientFilter::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os, indent);

  os << indent
     << "ResultArrayName:" << (this->ResultArrayName ? this->ResultArrayName : "Gradients") << endl;
  os << indent << "DivergenceArrayName:"
     << (this->DivergenceArrayName ? this->DivergenceArrayName : "Divergence") << endl;
  os << indent
     << "VorticityArrayName:" << (this->VorticityArrayName ? this->VorticityArrayName : "Vorticity")
     << endl;
  os << indent << "QCriterionArrayName:"
     << (this->QCriterionArrayName ? this->QCriterionArrayName : "Q-criterion") << endl;
  os << indent << "FasterApproximation:" << this->FasterApproximation << endl;
  os << indent << "ComputeGradient:" << this->ComputeGradient << endl;
  os << indent << "ComputeDivergence:" << this->ComputeDivergence << endl;
  os << indent << "ComputeVorticity:" << this->ComputeVorticity << endl;
  os << indent << "ComputeQCriterion:" << this->ComputeQCriterion << endl;
  os << indent << "ContributingCellOption:" << this->ContributingCellOption << endl;
  os << indent << "ReplacementValueOption:" << this->ReplacementValueOption << endl;
}

//------------------------------------------------------------------------------
void vtkGradientFilter::SetInputScalars(int fieldAssociation, const char* name)
{
  if ((fieldAssociation != vtkDataObject::FIELD_ASSOCIATION_POINTS) &&
    (fieldAssociation != vtkDataObject::FIELD_ASSOCIATION_CELLS) &&
    (fieldAssociation != vtkDataObject::FIELD_ASSOCIATION_POINTS_THEN_CELLS))
  {
    vtkErrorMacro("Input Array must be associated with points or cells.");
    return;
  }

  this->SetInputArrayToProcess(0, 0, 0, fieldAssociation, name);
}

//------------------------------------------------------------------------------
void vtkGradientFilter::SetInputScalars(int fieldAssociation, int fieldAttributeType)
{
  if ((fieldAssociation != vtkDataObject::FIELD_ASSOCIATION_POINTS) &&
    (fieldAssociation != vtkDataObject::FIELD_ASSOCIATION_CELLS) &&
    (fieldAssociation != vtkDataObject::FIELD_ASSOCIATION_POINTS_THEN_CELLS))
  {
    vtkErrorMacro("Input Array must be associated with points or cells.");
    return;
  }

  this->SetInputArrayToProcess(0, 0, 0, fieldAssociation, fieldAttributeType);
}

//------------------------------------------------------------------------------
int vtkGradientFilter::RequestUpdateExtent(vtkInformation* vtkNotUsed(request),
  vtkInformationVector** inputVector, vtkInformationVector* outputVector)
{
  // get the info objects
  vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
  vtkInformation* outInfo = outputVector->GetInformationObject(0);

  // Technically, this code is only correct for pieces extent types.  However,
  // since this class is pretty inefficient for data types that use 3D extents,
  // we'll punt on the ghost levels for them, too.
  int piece, numPieces, ghostLevels;

  piece = outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER());
  numPieces = outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES());
  ghostLevels = outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS());

  if (numPieces > 1)
  {
    ++ghostLevels;
  }

  inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER(), piece);
  inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES(), numPieces);
  inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS(), ghostLevels);
  inInfo->Set(vtkStreamingDemandDrivenPipeline::EXACT_EXTENT(), 1);

  return 1;
}

//------------------------------------------------------------------------------
int vtkGradientFilter::RequestData(vtkInformation* vtkNotUsed(request),
  vtkInformationVector** inputVector, vtkInformationVector* outputVector)
{
  vtkDebugMacro("RequestData");

  vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
  vtkInformation* outInfo = outputVector->GetInformationObject(0);

  vtkDataSet* input = vtkDataSet::SafeDownCast(inInfo->Get(vtkDataObject::DATA_OBJECT()));
  vtkDataSet* output = vtkDataSet::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));
  vtkDataArray* array = this->GetInputArrayToProcess(0, inputVector);

  if (input->GetNumberOfCells() == 0)
  {
    // need cells to compute the gradient so if we don't have cells. we can't compute anything.
    // if we have points and an array with values provide a warning letting the user know that
    // no gradient will be computed because of the lack of cells. otherwise the dataset is
    // assumed empty and we can skip providing a warning message to the user.
    if (input->GetNumberOfPoints() && array && array->GetNumberOfTuples())
    {
      vtkWarningMacro("Cannot compute gradient for datasets without cells");
    }
    output->ShallowCopy(input);
    return 1;
  }

  if (array == nullptr)
  {
    vtkErrorMacro("No input array. If this dataset is part of a composite dataset"
      << " check to make sure that all non-empty blocks have this array.");
    return 0;
  }
  if (array->GetNumberOfComponents() == 0)
  {
    vtkErrorMacro("Input array must have at least one component.");
    return 0;
  }

  // we can only compute vorticity and Q criterion if the input
  // array has 3 components. if we can't compute them because of
  // this we only mark internally the we aren't computing them
  // since we don't want to change the state of the filter.
  bool computeVorticity = this->ComputeVorticity != 0;
  bool computeDivergence = this->ComputeDivergence != 0;
  bool computeQCriterion = this->ComputeQCriterion != 0;
  if ((this->ComputeQCriterion || this->ComputeVorticity || this->ComputeDivergence) &&
    array->GetNumberOfComponents() != 3)
  {
    vtkWarningMacro("Input array must have exactly three components with "
      << "ComputeDivergence, ComputeVorticity or ComputeQCriterion flag enabled."
      << "Skipping divergence, vorticity and Q-criterion computation.");
    computeVorticity = false;
    computeQCriterion = false;
    computeDivergence = false;
  }

  int fieldAssociation;
  if (vtkGradientFilterHasArray(input->GetPointData(), array))
  {
    fieldAssociation = vtkDataObject::FIELD_ASSOCIATION_POINTS;
  }
  else if (vtkGradientFilterHasArray(input->GetCellData(), array))
  {
    fieldAssociation = vtkDataObject::FIELD_ASSOCIATION_CELLS;
  }
  else
  {
    vtkErrorMacro("Input arrays do not seem to be either point or cell arrays.");
    return 0;
  }

  output->CopyStructure(input);
  output->GetPointData()->PassData(input->GetPointData());
  output->GetCellData()->PassData(input->GetCellData());

  if (output->IsA("vtkImageData") || output->IsA("vtkStructuredGrid") ||
    output->IsA("vtkRectilinearGrid"))
  {
    this->ComputeRegularGridGradient(
      array, fieldAssociation, computeVorticity, computeQCriterion, computeDivergence, output);
  }
  else
  {
    this->ComputeUnstructuredGridGradient(array, fieldAssociation, input, computeVorticity,
      computeQCriterion, computeDivergence, output);
  }

  return 1;
}

//------------------------------------------------------------------------------
int vtkGradientFilter::ComputeUnstructuredGridGradient(vtkDataArray* array, int fieldAssociation,
  vtkDataSet* input, bool computeVorticity, bool computeQCriterion, bool computeDivergence,
  vtkDataSet* output)
{
  int arrayType = this->GetOutputArrayType(array);
  int numberOfInputComponents = array->GetNumberOfComponents();
  vtkSmartPointer<vtkDataArray> gradients = nullptr;
  if (this->ComputeGradient)
  {
    gradients.TakeReference(vtkDataArray::CreateDataArray(arrayType));
    gradients->SetNumberOfComponents(3 * numberOfInputComponents);
    gradients->SetNumberOfTuples(array->GetNumberOfTuples());
    switch (arrayType)
    {
      vtkFloatingPointTemplateMacro(
        Fill(gradients, static_cast<VTK_TT>(0), this->ReplacementValueOption));
    }
    if (this->ResultArrayName)
    {
      gradients->SetName(this->ResultArrayName);
    }
    else
    {
      gradients->SetName("Gradients");
    }
  }
  vtkSmartPointer<vtkDataArray> divergence = nullptr;
  if (computeDivergence)
  {
    divergence.TakeReference(vtkDataArray::CreateDataArray(arrayType));
    divergence->SetNumberOfTuples(array->GetNumberOfTuples());
    switch (arrayType)
    {
      vtkFloatingPointTemplateMacro(
        Fill(divergence, static_cast<VTK_TT>(0), this->ReplacementValueOption));
    }
    if (this->DivergenceArrayName)
    {
      divergence->SetName(this->DivergenceArrayName);
    }
    else
    {
      divergence->SetName("Divergence");
    }
  }
  vtkSmartPointer<vtkDataArray> vorticity;
  if (computeVorticity)
  {
    vorticity.TakeReference(vtkDataArray::CreateDataArray(arrayType));
    vorticity->SetNumberOfComponents(3);
    vorticity->SetNumberOfTuples(array->GetNumberOfTuples());
    switch (arrayType)
    {
      vtkFloatingPointTemplateMacro(
        Fill(vorticity, static_cast<VTK_TT>(0), this->ReplacementValueOption));
    }
    if (this->VorticityArrayName)
    {
      vorticity->SetName(this->VorticityArrayName);
    }
    else
    {
      vorticity->SetName("Vorticity");
    }
  }
  vtkSmartPointer<vtkDataArray> qCriterion;
  if (computeQCriterion)
  {
    qCriterion.TakeReference(vtkDataArray::CreateDataArray(arrayType));
    qCriterion->SetNumberOfTuples(array->GetNumberOfTuples());
    switch (arrayType)
    {
      vtkFloatingPointTemplateMacro(
        Fill(qCriterion, static_cast<VTK_TT>(0), this->ReplacementValueOption));
    }
    if (this->QCriterionArrayName)
    {
      qCriterion->SetName(this->QCriterionArrayName);
    }
    else
    {
      qCriterion->SetName("Q-criterion");
    }
  }

  int highestCellDimension = 0;
  if (this->ContributingCellOption == vtkGradientFilter::DataSetMax)
  {
    int maxDimension = input->IsA("vtkPolyData") == 1 ? 2 : 3;
    for (vtkIdType i = 0; i < input->GetNumberOfCells(); i++)
    {
      int dim = input->GetCell(i)->GetCellDimension();
      if (dim > highestCellDimension)
      {
        highestCellDimension = dim;
        if (highestCellDimension == maxDimension)
        {
          break;
        }
      }
    }
  }

  if (fieldAssociation == vtkDataObject::FIELD_ASSOCIATION_POINTS)
  {
    if (!this->FasterApproximation)
    {
      switch (arrayType)
      { // ok to use template macro here since we made the output arrays ourselves
        vtkFloatingPointTemplateMacro(ComputePointGradientsUG(input, array,
          (gradients == nullptr ? nullptr : static_cast<VTK_TT*>(gradients->GetVoidPointer(0))),
          numberOfInputComponents,
          (vorticity == nullptr ? nullptr : static_cast<VTK_TT*>(vorticity->GetVoidPointer(0))),
          (qCriterion == nullptr ? nullptr : static_cast<VTK_TT*>(qCriterion->GetVoidPointer(0))),
          (divergence == nullptr ? nullptr : static_cast<VTK_TT*>(divergence->GetVoidPointer(0))),
          highestCellDimension, this->ContributingCellOption));
      }
      if (gradients)
      {
        output->GetPointData()->AddArray(gradients);
      }
      if (divergence)
      {
        output->GetPointData()->AddArray(divergence);
      }
      if (vorticity)
      {
        output->GetPointData()->AddArray(vorticity);
      }
      if (qCriterion)
      {
        output->GetPointData()->AddArray(qCriterion);
      }
    }
    else // this->FasterApproximation
    {
      // The cell computation is faster and works off of point data anyway.  The
      // faster approximation is to use the cell algorithm and then convert the
      // result to point data.
      vtkSmartPointer<vtkDataArray> cellGradients = nullptr;
      if (gradients)
      {
        cellGradients.TakeReference(vtkDataArray::CreateDataArray(arrayType));
        cellGradients->SetName(gradients->GetName());
        cellGradients->SetNumberOfComponents(3 * array->GetNumberOfComponents());
        cellGradients->SetNumberOfTuples(input->GetNumberOfCells());
      }
      vtkSmartPointer<vtkDataArray> cellDivergence = nullptr;
      if (divergence)
      {
        cellDivergence.TakeReference(vtkDataArray::CreateDataArray(arrayType));
        cellDivergence->SetName(divergence->GetName());
        cellDivergence->SetNumberOfTuples(input->GetNumberOfCells());
      }
      vtkSmartPointer<vtkDataArray> cellVorticity = nullptr;
      if (vorticity)
      {
        cellVorticity.TakeReference(vtkDataArray::CreateDataArray(arrayType));
        cellVorticity->SetName(vorticity->GetName());
        cellVorticity->SetNumberOfComponents(3);
        cellVorticity->SetNumberOfTuples(input->GetNumberOfCells());
      }
      vtkSmartPointer<vtkDataArray> cellQCriterion = nullptr;
      if (qCriterion)
      {
        cellQCriterion.TakeReference(vtkDataArray::CreateDataArray(arrayType));
        cellQCriterion->SetName(qCriterion->GetName());
        cellQCriterion->SetNumberOfTuples(input->GetNumberOfCells());
      }

      switch (arrayType)
      { // ok to use template macro here since we made the output arrays ourselves
        vtkFloatingPointTemplateMacro(ComputeCellGradientsUG(input, array,
          (cellGradients == nullptr ? nullptr
                                    : static_cast<VTK_TT*>(cellGradients->GetVoidPointer(0))),
          numberOfInputComponents,
          (vorticity == nullptr ? nullptr : static_cast<VTK_TT*>(cellVorticity->GetVoidPointer(0))),
          (qCriterion == nullptr ? nullptr
                                 : static_cast<VTK_TT*>(cellQCriterion->GetVoidPointer(0))),
          (divergence == nullptr ? nullptr
                                 : static_cast<VTK_TT*>(cellDivergence->GetVoidPointer(0)))));
      }

      // We need to convert cell Array to points Array.
      vtkSmartPointer<vtkDataSet> dummy;
      dummy.TakeReference(input->NewInstance());
      dummy->CopyStructure(input);
      if (cellGradients)
      {
        dummy->GetCellData()->AddArray(cellGradients);
      }
      if (divergence)
      {
        dummy->GetCellData()->AddArray(cellDivergence);
      }
      if (vorticity)
      {
        dummy->GetCellData()->AddArray(cellVorticity);
      }
      if (qCriterion)
      {
        dummy->GetCellData()->AddArray(cellQCriterion);
      }

      vtkNew<vtkCellDataToPointData> cd2pd;
      cd2pd->SetInputData(dummy);
      cd2pd->PassCellDataOff();
      cd2pd->SetContributingCellOption(this->ContributingCellOption);
      cd2pd->Update();

      // Set the gradients array in the output and cleanup.
      if (gradients)
      {
        output->GetPointData()->AddArray(
          cd2pd->GetOutput()->GetPointData()->GetArray(gradients->GetName()));
      }
      if (qCriterion)
      {
        output->GetPointData()->AddArray(
          cd2pd->GetOutput()->GetPointData()->GetArray(qCriterion->GetName()));
      }
      if (divergence)
      {
        output->GetPointData()->AddArray(
          cd2pd->GetOutput()->GetPointData()->GetArray(divergence->GetName()));
      }
      if (vorticity)
      {
        output->GetPointData()->AddArray(
          cd2pd->GetOutput()->GetPointData()->GetArray(vorticity->GetName()));
      }
    }
  }
  else // fieldAssociation == vtkDataObject::FIELD_ASSOCIATION_CELLS
  {
    // We need to convert cell Array to points Array.
    vtkSmartPointer<vtkDataSet> dummy;
    dummy.TakeReference(input->NewInstance());
    dummy->CopyStructure(input);
    dummy->GetCellData()->SetScalars(array);

    vtkNew<vtkCellDataToPointData> cd2pd;
    cd2pd->SetInputData(dummy);
    cd2pd->PassCellDataOff();
    cd2pd->SetContributingCellOption(this->ContributingCellOption);
    cd2pd->Update();
    vtkDataArray* pointScalars = cd2pd->GetOutput()->GetPointData()->GetScalars();
    pointScalars->Register(this);

    switch (arrayType)
    { // ok to use template macro here since we made the output arrays ourselves
      vtkFloatingPointTemplateMacro(ComputeCellGradientsUG(input, pointScalars,
        (gradients == nullptr ? nullptr : static_cast<VTK_TT*>(gradients->GetVoidPointer(0))),
        numberOfInputComponents,
        (vorticity == nullptr ? nullptr : static_cast<VTK_TT*>(vorticity->GetVoidPointer(0))),
        (qCriterion == nullptr ? nullptr : static_cast<VTK_TT*>(qCriterion->GetVoidPointer(0))),
        (divergence == nullptr ? nullptr : static_cast<VTK_TT*>(divergence->GetVoidPointer(0)))));
    }

    if (gradients)
    {
      output->GetCellData()->AddArray(gradients);
    }
    if (vorticity)
    {
      output->GetCellData()->AddArray(vorticity);
    }
    if (divergence)
    {
      output->GetCellData()->AddArray(divergence);
    }
    if (qCriterion)
    {
      output->GetCellData()->AddArray(qCriterion);
    }
    pointScalars->UnRegister(this);
  }

  return 1;
}

//------------------------------------------------------------------------------
int vtkGradientFilter::ComputeRegularGridGradient(vtkDataArray* array, int fieldAssociation,
  bool computeVorticity, bool computeQCriterion, bool computeDivergence, vtkDataSet* output)
{
  int arrayType = this->GetOutputArrayType(array);
  int numberOfInputComponents = array->GetNumberOfComponents();
  vtkSmartPointer<vtkDataArray> gradients = nullptr;
  if (this->ComputeGradient)
  {
    gradients.TakeReference(vtkDataArray::CreateDataArray(arrayType));
    gradients->SetNumberOfComponents(3 * numberOfInputComponents);
    gradients->SetNumberOfTuples(array->GetNumberOfTuples());
    if (this->ResultArrayName)
    {
      gradients->SetName(this->ResultArrayName);
    }
    else
    {
      gradients->SetName("Gradients");
    }
  }
  vtkSmartPointer<vtkDataArray> divergence = nullptr;
  if (computeDivergence)
  {
    divergence.TakeReference(vtkDataArray::CreateDataArray(arrayType));
    divergence->SetNumberOfTuples(array->GetNumberOfTuples());
    if (this->DivergenceArrayName)
    {
      divergence->SetName(this->DivergenceArrayName);
    }
    else
    {
      divergence->SetName("Divergence");
    }
  }
  vtkSmartPointer<vtkDataArray> vorticity;
  if (computeVorticity)
  {
    vorticity.TakeReference(vtkDataArray::CreateDataArray(arrayType));
    vorticity->SetNumberOfComponents(3);
    vorticity->SetNumberOfTuples(array->GetNumberOfTuples());
    if (this->VorticityArrayName)
    {
      vorticity->SetName(this->VorticityArrayName);
    }
    else
    {
      vorticity->SetName("Vorticity");
    }
  }
  vtkSmartPointer<vtkDataArray> qCriterion;
  if (computeQCriterion)
  {
    qCriterion.TakeReference(vtkDataArray::CreateDataArray(arrayType));
    qCriterion->SetNumberOfTuples(array->GetNumberOfTuples());
    if (this->QCriterionArrayName)
    {
      qCriterion->SetName(this->QCriterionArrayName);
    }
    else
    {
      qCriterion->SetName("Q-criterion");
    }
  }

  if (vtkStructuredGrid* structuredGrid = vtkStructuredGrid::SafeDownCast(output))
  {
    switch (arrayType)
    { // ok to use template macro here since we made the output arrays ourselves
      vtkFloatingPointTemplateMacro(ComputeGradientsSG(structuredGrid, array,
        (gradients == nullptr ? nullptr : static_cast<VTK_TT*>(gradients->GetVoidPointer(0))),
        numberOfInputComponents, fieldAssociation,
        (vorticity == nullptr ? nullptr : static_cast<VTK_TT*>(vorticity->GetVoidPointer(0))),
        (qCriterion == nullptr ? nullptr : static_cast<VTK_TT*>(qCriterion->GetVoidPointer(0))),
        (divergence == nullptr ? nullptr : static_cast<VTK_TT*>(divergence->GetVoidPointer(0)))));
    }
  }
  else if (vtkImageData* imageData = vtkImageData::SafeDownCast(output))
  {
    switch (arrayType)
    { // ok to use template macro here since we made the output arrays ourselves
      vtkFloatingPointTemplateMacro(ComputeGradientsSG(imageData, array,
        (gradients == nullptr ? nullptr : static_cast<VTK_TT*>(gradients->GetVoidPointer(0))),
        numberOfInputComponents, fieldAssociation,
        (vorticity == nullptr ? nullptr : static_cast<VTK_TT*>(vorticity->GetVoidPointer(0))),
        (qCriterion == nullptr ? nullptr : static_cast<VTK_TT*>(qCriterion->GetVoidPointer(0))),
        (divergence == nullptr ? nullptr : static_cast<VTK_TT*>(divergence->GetVoidPointer(0)))));
    }
  }
  else if (vtkRectilinearGrid* rectilinearGrid = vtkRectilinearGrid::SafeDownCast(output))
  {
    switch (arrayType)
    { // ok to use template macro here since we made the output arrays ourselves
      vtkFloatingPointTemplateMacro(ComputeGradientsSG(rectilinearGrid, array,
        (gradients == nullptr ? nullptr : static_cast<VTK_TT*>(gradients->GetVoidPointer(0))),
        numberOfInputComponents, fieldAssociation,
        (vorticity == nullptr ? nullptr : static_cast<VTK_TT*>(vorticity->GetVoidPointer(0))),
        (qCriterion == nullptr ? nullptr : static_cast<VTK_TT*>(qCriterion->GetVoidPointer(0))),
        (divergence == nullptr ? nullptr : static_cast<VTK_TT*>(divergence->GetVoidPointer(0)))));
    }
  }
  if (fieldAssociation == vtkDataObject::FIELD_ASSOCIATION_POINTS)
  {
    if (gradients)
    {
      output->GetPointData()->AddArray(gradients);
    }
    if (vorticity)
    {
      output->GetPointData()->AddArray(vorticity);
    }
    if (qCriterion)
    {
      output->GetPointData()->AddArray(qCriterion);
    }
    if (divergence)
    {
      output->GetPointData()->AddArray(divergence);
    }
  }
  else if (fieldAssociation == vtkDataObject::FIELD_ASSOCIATION_CELLS)
  {
    if (gradients)
    {
      output->GetCellData()->AddArray(gradients);
    }
    if (vorticity)
    {
      output->GetCellData()->AddArray(vorticity);
    }
    if (qCriterion)
    {
      output->GetCellData()->AddArray(qCriterion);
    }
    if (divergence)
    {
      output->GetCellData()->AddArray(divergence);
    }
  }
  else
  {
    vtkErrorMacro("Bad fieldAssociation value " << fieldAssociation << endl);
  }

  return 1;
}

//------------------------------------------------------------------------------
int vtkGradientFilter::GetOutputArrayType(vtkDataArray* array)
{
  int retType = VTK_DOUBLE;
  switch (array->GetDataType())
  {
    vtkTemplateMacro(retType = GetOutputDataType(static_cast<VTK_TT>(0)));
  }
  return retType;
}

namespace
{
//------------------------------------------------------------------------------
template <class data_type>
void ComputePointGradientsUG(vtkDataSet* structure, vtkDataArray* array, data_type* gradients,
  int numberOfInputComponents, data_type* vorticity, data_type* qCriterion, data_type* divergence,
  int highestCellDimension, int contributingCellOption)
{
  vtkNew<vtkIdList> currentPoint;
  currentPoint->SetNumberOfIds(1);
  vtkNew<vtkIdList> cellsOnPoint;

  vtkIdType numpts = structure->GetNumberOfPoints();

  int numberOfOutputComponents = 3 * numberOfInputComponents;
  std::vector<data_type> g(numberOfOutputComponents);

  // if we are doing patches for contributing cell dimensions we want to keep track of
  // the maximum expected dimension so we can exit out of the check loop quicker
  const int maxCellDimension = structure->IsA("vtkPolyData") ? 2 : 3;

  for (vtkIdType point = 0; point < numpts; point++)
  {
    currentPoint->SetId(0, point);
    double pointcoords[3];
    structure->GetPoint(point, pointcoords);
    // Get all cells touching this point.
    structure->GetCellNeighbors(-1, currentPoint, cellsOnPoint);
    vtkIdType numCellNeighbors = cellsOnPoint->GetNumberOfIds();

    for (int i = 0; i < numberOfOutputComponents; i++)
    {
      g[i] = 0;
    }

    if (contributingCellOption == vtkGradientFilter::Patch)
    {
      highestCellDimension = 0;
      for (vtkIdType neighbor = 0; neighbor < numCellNeighbors; neighbor++)
      {
        int cellDimension = structure->GetCell(cellsOnPoint->GetId(neighbor))->GetCellDimension();
        if (cellDimension > highestCellDimension)
        {
          highestCellDimension = cellDimension;
          if (highestCellDimension == maxCellDimension)
          {
            break;
          }
        }
      }
    }
    vtkIdType numValidCellNeighbors = 0;

    // Iterate on all cells and find all points connected to current point
    // by an edge.
    for (vtkIdType neighbor = 0; neighbor < numCellNeighbors; neighbor++)
    {
      vtkCell* cell = structure->GetCell(cellsOnPoint->GetId(neighbor));
      if (cell->GetCellDimension() >= highestCellDimension)
      {
        int subId;
        double parametricCoord[3];
        if (GetCellParametricData(point, pointcoords, cell, subId, parametricCoord))
        {
          numValidCellNeighbors++;
          for (int inputComponent = 0; inputComponent < numberOfInputComponents; inputComponent++)
          {
            int numberOfCellPoints = cell->GetNumberOfPoints();
            std::vector<double> values(numberOfCellPoints);
            // Get values of Array at cell points.
            for (int i = 0; i < numberOfCellPoints; i++)
            {
              values[i] = array->GetComponent(cell->GetPointId(i), inputComponent);
            }

            double derivative[3];
            // Get derivative of cell at point.
            cell->Derivatives(subId, parametricCoord, &values[0], 1, derivative);

            g[inputComponent * 3] += static_cast<data_type>(derivative[0]);
            g[inputComponent * 3 + 1] += static_cast<data_type>(derivative[1]);
            g[inputComponent * 3 + 2] += static_cast<data_type>(derivative[2]);
          } // iterating over Components
        }   // if(GetCellParametricData())
      }     // if(cell->GetCellDimension () >= highestCellDimension
    }       // iterating over neighbors

    if (numValidCellNeighbors > 0)
    {
      for (int i = 0; i < 3 * numberOfInputComponents; i++)
      {
        g[i] /= numValidCellNeighbors;
      }

      if (vorticity)
      {
        ComputeVorticityFromGradient(&g[0], vorticity + 3 * point);
      }
      if (qCriterion)
      {
        ComputeQCriterionFromGradient(&g[0], qCriterion + point);
      }
      if (divergence)
      {
        ComputeDivergenceFromGradient(&g[0], divergence + point);
      }
      if (gradients)
      {
        for (int i = 0; i < numberOfOutputComponents; i++)
        {
          gradients[point * numberOfOutputComponents + i] = g[i];
        }
      }
    }
  } // iterating over points in grid
}

//------------------------------------------------------------------------------
int GetCellParametricData(
  vtkIdType pointId, double pointCoord[3], vtkCell* cell, int& subId, double parametricCoord[3])
{
  // Watch out for degenerate cells.  They make the derivative calculation
  // fail.
  vtkIdList* pointIds = cell->GetPointIds();
  int timesPointRegistered = 0;
  for (int i = 0; i < pointIds->GetNumberOfIds(); i++)
  {
    if (pointId == pointIds->GetId(i))
    {
      timesPointRegistered++;
    }
  }
  if (timesPointRegistered != 1)
  {
    // The cell should have the point exactly once.  Not good.
    return 0;
  }

  double dummy;
  int numpoints = cell->GetNumberOfPoints();
  std::vector<double> values(numpoints);
  // Get parametric position of point.
  cell->EvaluatePosition(
    pointCoord, nullptr, subId, parametricCoord, dummy, &values[0] /*Really another dummy.*/);

  return 1;
}

//------------------------------------------------------------------------------
template <class data_type>
void ComputeCellGradientsUG(vtkDataSet* structure, vtkDataArray* array, data_type* gradients,
  int numberOfInputComponents, data_type* vorticity, data_type* qCriterion, data_type* divergence)
{
  vtkIdType numcells = structure->GetNumberOfCells();
  std::vector<double> values(8);
  std::vector<data_type> cellGradients(3 * numberOfInputComponents);
  for (vtkIdType cellid = 0; cellid < numcells; cellid++)
  {
    vtkCell* cell = structure->GetCell(cellid);
    int subId;
    double cellCenter[3];
    subId = cell->GetParametricCenter(cellCenter);

    int numpoints = cell->GetNumberOfPoints();
    if (static_cast<size_t>(numpoints) > values.size())
    {
      values.resize(numpoints);
    }
    double derivative[3];
    for (int inputComponent = 0; inputComponent < numberOfInputComponents; inputComponent++)
    {
      for (int i = 0; i < numpoints; i++)
      {
        values[i] = array->GetComponent(cell->GetPointId(i), inputComponent);
      }

      cell->Derivatives(subId, cellCenter, &values[0], 1, derivative);
      cellGradients[inputComponent * 3] = static_cast<data_type>(derivative[0]);
      cellGradients[inputComponent * 3 + 1] = static_cast<data_type>(derivative[1]);
      cellGradients[inputComponent * 3 + 2] = static_cast<data_type>(derivative[2]);
    }
    if (gradients)
    {
      for (int i = 0; i < 3 * numberOfInputComponents; i++)
      {
        gradients[cellid * 3 * numberOfInputComponents + i] = cellGradients[i];
      }
    }
    if (vorticity)
    {
      ComputeVorticityFromGradient(&cellGradients[0], vorticity + 3 * cellid);
    }
    if (qCriterion)
    {
      ComputeQCriterionFromGradient(&cellGradients[0], qCriterion + cellid);
    }
    if (divergence)
    {
      ComputeDivergenceFromGradient(&cellGradients[0], divergence + cellid);
    }
  }
}

//------------------------------------------------------------------------------
template <class Grid, class data_type>
void ComputeGradientsSG(Grid output, vtkDataArray* array, data_type* gradients,
  int numberOfInputComponents, int fieldAssociation, data_type* vorticity, data_type* qCriterion,
  data_type* divergence)
{
  int idx, idx2, inputComponent;
  double xp[3], xm[3], factor;
  xp[0] = xp[1] = xp[2] = xm[0] = xm[1] = xm[2] = factor = 0;
  double xxi, yxi, zxi, xeta, yeta, zeta, xzeta, yzeta, zzeta;
  yxi = zxi = xeta = yeta = zeta = xzeta = yzeta = zzeta = 0;
  double aj, xix, xiy, xiz, etax, etay, etaz, zetax, zetay, zetaz;
  xix = xiy = xiz = etax = etay = etaz = zetax = zetay = zetaz = 0;
  // for finite differencing -- the values on the "plus" side and
  // "minus" side of the point to be computed at
  std::vector<double> plusvalues(numberOfInputComponents);
  std::vector<double> minusvalues(numberOfInputComponents);

  std::vector<double> dValuesdXi(numberOfInputComponents);
  std::vector<double> dValuesdEta(numberOfInputComponents);
  std::vector<double> dValuesdZeta(numberOfInputComponents);
  std::vector<data_type> localGradients(numberOfInputComponents * 3);

  int dims[3];
  output->GetDimensions(dims);
  if (fieldAssociation == vtkDataObject::FIELD_ASSOCIATION_CELLS)
  {
    // reduce the dimensions by 1 for cells
    for (int i = 0; i < 3; i++)
    {
      dims[i]--;
    }
  }
  int ijsize = dims[0] * dims[1];

  for (int k = 0; k < dims[2]; k++)
  {
    for (int j = 0; j < dims[1]; j++)
    {
      for (int i = 0; i < dims[0]; i++)
      {
        //  Xi derivatives.
        if (dims[0] == 1) // 2D in this direction
        {
          factor = 1.0;
          for (int ii = 0; ii < 3; ii++)
          {
            xp[ii] = xm[ii] = 0.0;
          }
          xp[0] = 1.0;
          for (inputComponent = 0; inputComponent < numberOfInputComponents; inputComponent++)
          {
            plusvalues[inputComponent] = minusvalues[inputComponent] = 0;
          }
        }
        else if (i == 0)
        {
          factor = 1.0;
          idx = (i + 1) + j * dims[0] + k * ijsize;
          idx2 = i + j * dims[0] + k * ijsize;
          GetGridEntityCoordinate(output, fieldAssociation, idx, xp);
          GetGridEntityCoordinate(output, fieldAssociation, idx2, xm);
          for (inputComponent = 0; inputComponent < numberOfInputComponents; inputComponent++)
          {
            plusvalues[inputComponent] = array->GetComponent(idx, inputComponent);
            minusvalues[inputComponent] = array->GetComponent(idx2, inputComponent);
          }
        }
        else if (i == (dims[0] - 1))
        {
          factor = 1.0;
          idx = i + j * dims[0] + k * ijsize;
          idx2 = i - 1 + j * dims[0] + k * ijsize;
          GetGridEntityCoordinate(output, fieldAssociation, idx, xp);
          GetGridEntityCoordinate(output, fieldAssociation, idx2, xm);
          for (inputComponent = 0; inputComponent < numberOfInputComponents; inputComponent++)
          {
            plusvalues[inputComponent] = array->GetComponent(idx, inputComponent);
            minusvalues[inputComponent] = array->GetComponent(idx2, inputComponent);
          }
        }
        else
        {
          factor = 0.5;
          idx = (i + 1) + j * dims[0] + k * ijsize;
          idx2 = (i - 1) + j * dims[0] + k * ijsize;
          GetGridEntityCoordinate(output, fieldAssociation, idx, xp);
          GetGridEntityCoordinate(output, fieldAssociation, idx2, xm);
          for (inputComponent = 0; inputComponent < numberOfInputComponents; inputComponent++)
          {
            plusvalues[inputComponent] = array->GetComponent(idx, inputComponent);
            minusvalues[inputComponent] = array->GetComponent(idx2, inputComponent);
          }
        }

        xxi = factor * (xp[0] - xm[0]);
        yxi = factor * (xp[1] - xm[1]);
        zxi = factor * (xp[2] - xm[2]);
        for (inputComponent = 0; inputComponent < numberOfInputComponents; inputComponent++)
        {
          dValuesdXi[inputComponent] =
            factor * (plusvalues[inputComponent] - minusvalues[inputComponent]);
        }

        //  Eta derivatives.
        if (dims[1] == 1) // 2D in this direction
        {
          factor = 1.0;
          for (int ii = 0; ii < 3; ii++)
          {
            xp[ii] = xm[ii] = 0.0;
          }
          xp[1] = 1.0;
          for (inputComponent = 0; inputComponent < numberOfInputComponents; inputComponent++)
          {
            plusvalues[inputComponent] = minusvalues[inputComponent] = 0;
          }
        }
        else if (j == 0)
        {
          factor = 1.0;
          idx = i + (j + 1) * dims[0] + k * ijsize;
          idx2 = i + j * dims[0] + k * ijsize;
          GetGridEntityCoordinate(output, fieldAssociation, idx, xp);
          GetGridEntityCoordinate(output, fieldAssociation, idx2, xm);
          for (inputComponent = 0; inputComponent < numberOfInputComponents; inputComponent++)
          {
            plusvalues[inputComponent] = array->GetComponent(idx, inputComponent);
            minusvalues[inputComponent] = array->GetComponent(idx2, inputComponent);
          }
        }
        else if (j == (dims[1] - 1))
        {
          factor = 1.0;
          idx = i + j * dims[0] + k * ijsize;
          idx2 = i + (j - 1) * dims[0] + k * ijsize;
          GetGridEntityCoordinate(output, fieldAssociation, idx, xp);
          GetGridEntityCoordinate(output, fieldAssociation, idx2, xm);
          for (inputComponent = 0; inputComponent < numberOfInputComponents; inputComponent++)
          {
            plusvalues[inputComponent] = array->GetComponent(idx, inputComponent);
            minusvalues[inputComponent] = array->GetComponent(idx2, inputComponent);
          }
        }
        else
        {
          factor = 0.5;
          idx = i + (j + 1) * dims[0] + k * ijsize;
          idx2 = i + (j - 1) * dims[0] + k * ijsize;
          GetGridEntityCoordinate(output, fieldAssociation, idx, xp);
          GetGridEntityCoordinate(output, fieldAssociation, idx2, xm);
          for (inputComponent = 0; inputComponent < numberOfInputComponents; inputComponent++)
          {
            plusvalues[inputComponent] = array->GetComponent(idx, inputComponent);
            minusvalues[inputComponent] = array->GetComponent(idx2, inputComponent);
          }
        }

        xeta = factor * (xp[0] - xm[0]);
        yeta = factor * (xp[1] - xm[1]);
        zeta = factor * (xp[2] - xm[2]);
        for (inputComponent = 0; inputComponent < numberOfInputComponents; inputComponent++)
        {
          dValuesdEta[inputComponent] =
            factor * (plusvalues[inputComponent] - minusvalues[inputComponent]);
        }

        //  Zeta derivatives.
        if (dims[2] == 1) // 2D in this direction
        {
          factor = 1.0;
          for (int ii = 0; ii < 3; ii++)
          {
            xp[ii] = xm[ii] = 0.0;
          }
          for (inputComponent = 0; inputComponent < numberOfInputComponents; inputComponent++)
          {
            plusvalues[inputComponent] = minusvalues[inputComponent] = 0;
          }
          xp[2] = 1.0;
        }
        else if (k == 0)
        {
          factor = 1.0;
          idx = i + j * dims[0] + (k + 1) * ijsize;
          idx2 = i + j * dims[0] + k * ijsize;
          GetGridEntityCoordinate(output, fieldAssociation, idx, xp);
          GetGridEntityCoordinate(output, fieldAssociation, idx2, xm);
          for (inputComponent = 0; inputComponent < numberOfInputComponents; inputComponent++)
          {
            plusvalues[inputComponent] = array->GetComponent(idx, inputComponent);
            minusvalues[inputComponent] = array->GetComponent(idx2, inputComponent);
          }
        }
        else if (k == (dims[2] - 1))
        {
          factor = 1.0;
          idx = i + j * dims[0] + k * ijsize;
          idx2 = i + j * dims[0] + (k - 1) * ijsize;
          GetGridEntityCoordinate(output, fieldAssociation, idx, xp);
          GetGridEntityCoordinate(output, fieldAssociation, idx2, xm);
          for (inputComponent = 0; inputComponent < numberOfInputComponents; inputComponent++)
          {
            plusvalues[inputComponent] = array->GetComponent(idx, inputComponent);
            minusvalues[inputComponent] = array->GetComponent(idx2, inputComponent);
          }
        }
        else
        {
          factor = 0.5;
          idx = i + j * dims[0] + (k + 1) * ijsize;
          idx2 = i + j * dims[0] + (k - 1) * ijsize;
          GetGridEntityCoordinate(output, fieldAssociation, idx, xp);
          GetGridEntityCoordinate(output, fieldAssociation, idx2, xm);
          for (inputComponent = 0; inputComponent < numberOfInputComponents; inputComponent++)
          {
            plusvalues[inputComponent] = array->GetComponent(idx, inputComponent);
            minusvalues[inputComponent] = array->GetComponent(idx2, inputComponent);
          }
        }

        xzeta = factor * (xp[0] - xm[0]);
        yzeta = factor * (xp[1] - xm[1]);
        zzeta = factor * (xp[2] - xm[2]);
        for (inputComponent = 0; inputComponent < numberOfInputComponents; inputComponent++)
        {
          dValuesdZeta[inputComponent] =
            factor * (plusvalues[inputComponent] - minusvalues[inputComponent]);
        }

        // Now calculate the Jacobian.  Grids occasionally have
        // singularities, or points where the Jacobian is infinite (the
        // inverse is zero).  For these cases, we'll set the Jacobian to
        // zero, which will result in a zero derivative.
        //
        aj = xxi * yeta * zzeta + yxi * zeta * xzeta + zxi * xeta * yzeta - zxi * yeta * xzeta -
          yxi * xeta * zzeta - xxi * zeta * yzeta;
        if (aj != 0.0)
        {
          aj = 1. / aj;
        }

        //  Xi metrics.
        xix = aj * (yeta * zzeta - zeta * yzeta);
        xiy = -aj * (xeta * zzeta - zeta * xzeta);
        xiz = aj * (xeta * yzeta - yeta * xzeta);

        //  Eta metrics.
        etax = -aj * (yxi * zzeta - zxi * yzeta);
        etay = aj * (xxi * zzeta - zxi * xzeta);
        etaz = -aj * (xxi * yzeta - yxi * xzeta);

        //  Zeta metrics.
        zetax = aj * (yxi * zeta - zxi * yeta);
        zetay = -aj * (xxi * zeta - zxi * xeta);
        zetaz = aj * (xxi * yeta - yxi * xeta);

        // Finally compute the actual derivatives
        idx = i + j * dims[0] + k * ijsize;
        for (inputComponent = 0; inputComponent < numberOfInputComponents; inputComponent++)
        {
          localGradients[inputComponent * 3] =
            static_cast<data_type>(xix * dValuesdXi[inputComponent] +
              etax * dValuesdEta[inputComponent] + zetax * dValuesdZeta[inputComponent]);

          localGradients[inputComponent * 3 + 1] =
            static_cast<data_type>(xiy * dValuesdXi[inputComponent] +
              etay * dValuesdEta[inputComponent] + zetay * dValuesdZeta[inputComponent]);

          localGradients[inputComponent * 3 + 2] =
            static_cast<data_type>(xiz * dValuesdXi[inputComponent] +
              etaz * dValuesdEta[inputComponent] + zetaz * dValuesdZeta[inputComponent]);
        }

        if (gradients)
        {
          for (int ii = 0; ii < 3 * numberOfInputComponents; ii++)
          {
            gradients[idx * numberOfInputComponents * 3 + ii] = localGradients[ii];
          }
        }
        if (vorticity)
        {
          ComputeVorticityFromGradient(&localGradients[0], vorticity + 3 * idx);
        }
        if (qCriterion)
        {
          ComputeQCriterionFromGradient(&localGradients[0], qCriterion + idx);
        }
        if (divergence)
        {
          ComputeDivergenceFromGradient(&localGradients[0], divergence + idx);
        }
      }
    }
  }
}

} // end anonymous namespace
